Mechanisms of budding of nanoscale particles 

through lipid bilayers 

Teresa Ruiz-Herrero,*'^ Enrique Velasco,^ and Michael F. Hagan*'* 

Departamento de Fisica Teorica de la Materia Condensada, Universidad Autdnoma de Madrid, 
Madrid, and Department of Physics, Brandeis University, Waltham, MA 

^ E-mail: _teresa.ruiz@uam.es ; hagan@brandeis.edu] 

O 

(N 

X) 
D 

tin 

u 

d 
• ^ 

cr 



O 



X 



*To whom correspondence should be addressed 
^Universidad Autonoma de Madrid 
^Brandeis University 



1 



Abstract 



We examine the budding of a nanoscale particle through a lipid bilayer using molecular dynamics 
simulations, free energy calculations, and an elastic theory, with the aim of determining the extent 
to which equilibrium elasticity theory can describe the factors that control the mechanism and effi- 
ciency of budding. The particle is a smooth sphere which experiences attractive interactions to the 
lipid head groups. Depending on the parameters, we observe four classes of dynamical trajectories: 
particle adhesion to the membrane, stalled partially wrapped states, budding followed by scission, and 
membrane rupture. In most regions of parameter space we find that the elastic theory agrees nearly 
quantitatively with the simulated phase behavior as a function of adhesion strength, membrane bending 
rigidity, and particle radius. However, at parameter values near the transition between particle adhesion 
and budding, we observe long-lived partially wrapped states which are not captured by existing elastic 
theories. These states could constrain the accessible system parameters for those enveloped viruses or 
drug delivery vehicles which rely on exo- or endocytosis for membrane transport. 



Introduction 

The mechanisms by which nanoscale particles 
cross cell membranes and the factors that con- 
trol their uptake are essential questions for cellular 
physiology and modern biomedicine. Regulating 
the uptake (endocytosis) of nanoparticles is impor- 
tant for nanomedicine applications and for predict- 
ing nanoparticle toxicity .l^^ Similarly, during the 
replication of many viruses an assembled nucleo- 
capsid buds through the cell membrane, simulta- 
neously exiting the cell and acquiring a membrane 
coating of host origin. Although endocytosis,^ 
viral budding,!^ and scission of budded viruses^ 
can be actively driven or assisted by cell machin- 
ery, both nanoparticle uptake and at least some as- 
pects of budding of viruses or viral proteins'^ can 
occur passively(without cell machinery or ATP 
hydrolysis).Q2HiII Furthermore, evidence suggests 
that some viruses do o r can undergo passive bud- 
ding in vivo (e.gJ^^^^. It is therefore important 
to establish the aspects of particle budding which 
are generic to passive transport and thus under- 
lie all forms of particle uptake or egress. In this 
paper we use elastic theory, moleculary dynam- 
ics (MD) simulations, and free energy calculations 
to characterize the dynamics and thermodynam- 
ics of the process by which a particle adheres to 
a membrane, is passively engulfed, and then spon- 
taneously separates. 

Previous works first studied the equilibrium con- 
figurations of budding through a vesicle or infi- 
nite membrane as a function of membrane rigidity, 
particle size, and membrane-particle adhesion en- 



ergy using elasticity theory-^^^^^^ Subsequent stud- 
ies used simulations to address further aspects 
of the problem, including Monte Carlo simula- 
tions on a randomly triangulated surface repre- 
sentation of a vesicle and molecular dynamics 
(MD) on a coarse-grained lipid model to in- 
vestigate wrapping of charged particles, density 
functional theory to study the relationship between 
particle hydrophobicity and wrapping, and dis- 
sipative particle dynamics (DPD) to study wrap- 
ping of a particle by a inhomogeneous bilayer"^ 
and wrapping behavior of ligand-coated nanopar- 
ticles.^^ Recently, the wrapping behavior of ellip- 
soidal particles has been studied via DPD^^ and 
MD simulations 

While all of these treatments show that the adhe- 
sion energy required for wrapping depends on par- 
ticle properties and membrane composition, there 
has not been a thorough comparison of predictions 
of elasticity theory with the results of more sophis- 
ticated computational models. In this work our 
primary objective is to understand the extent to 
which simplified elastic models can describe the 
thermodynamics and/or dynamics of particle up- 
take. To focus on aspects generic to all forms of 
exo- or endocytosis, we consider a minimal model 
in which the membrane is treated as a bilayer of 
homogeneous composition and the particle is pre- 
assembled (as in th e cas e of nanoparticles or, e.g. 
type-D retroviruses'^^'^, and is spherically sym- 
metric. Thus, in this work we do not consider 
the effects of membrane inhomogeneity (i.e. lipid 
rafts) or the association of viral membrane pro- 
teins.-^— We compare the predictions of the elastic 
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model to results of dynamical simulations and 
free energy calculations. We find that the phase 
behavior predicted by the two descriptions agrees 
nearly quantitatively in most regions of parame- 
ter space, but there are important dynamical dif- 
ferences at parameter values near the transition be- 
tween no uptake and particle budding. In particu- 
lar, we identify a partially wrapped state which we 
show to be metastable. 

Methods 

The Membrane Model 

We model the amphiphilic lipids comprising the 
membrane with a coarse grained implicit solvent 
model from Cooke et in which each am- 
phiphile is represented by one head bead and two 
tail beads that interact via repulsive WCA poten- 
tials, ^9 Eq.(l) 
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with rc = and b is chosen to ensure an effec- 
tive cylindrical lipid shape: Z7head-head = ^head-tail = 
0.95(7 and 

^head-tail — where o wiU turn out 
to be the typical distance between beads within a 
model lipid molecule. 

The beads belonging to a given lipid are con- 
nected through FENE bonds (Eq.(2))'^ and the 
linearity of the molecule is achieved via a har- 
monic spring with rest length 4o between the first 
and the third bead, Eq.(3) 
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Since this is an implicit solvent model, the hy- 
drophobicity is represented by an attractive inter- 
action, Eq.(4), between all tail beads. 
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Figure 1 : The properties of the membrane are easily 
tuned by (Oc- (a) Bending rigidity K (in log scale) and 
(b) areal density T] of lipids as functions of (Oc- The 
values in (a) are from Ref . and the values in (b) were 
calculated from our simulations. 

This model allows the formation of bilayers with 
physical properties such as fluidity, area per 
molecule and bending rigidity that are easily tuned 
via (Oc- Moreover, diffusivity within the mem- 
brane, density, and bending rigidity are in good 
agreement with values of theseparameters mea- 
sured for biological membranes'^ (Figure 1) 

Membrane-particle interaction 

As noted in the introduction, the systems we 
have in mind include synthetic nanoparticles or 
viral particles which bud through attractive in- 
teractions with lipid membranes. These interac- 
tions can arise in part from electrostatic inter- 
actions between charged lipid head groups and 
charges on the nanoparticle surface or capsid ex- 
terior (e.g. basic residues on the matrix protein 
in retroviruses'^. A second source of interac- 
tion can be protein mediated, including binding of 
nanoparticle-functionalized ligands to membrane 
receptors or insertion of hydrophobic tails on cap- 
sid proteins into the membrane.^ ^Finally, trans- 
membrane viral 'spike' proteins can drive or facil- 
itate budding. Importantly, each of these forms of 
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interactions is short ranged. Receptor-ligand and 
spike protein-virus interactions operate on length 
scales of A to nm; similarly, at physiological con- 
ditions of 100 mM salt electrostatic interactions 
have a Debye screening length of 1 nm. Thus, 
to keep our analysis general, we consider a short 
range attractive interaction between our model 
particle and head groups. In particular, we rep- 
resent the combination of excluded volume and at- 
tractive interactions between the particle and head 
groups with a shifted Lennard Jones potential: 
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with £ a free parameter that controls the 
membrane-particle interaction strength, s = R — 

(7/2, the particle radius, Vcut = 4£ (^^)^^- (^) 

rph = 3.5a, and e - Vcut the depth of the attractive in- 
teraction between the particle and the membrane. 

The particle experiences only excluded volume 
interactions with the tail groups, which are mod- 
eled with a shifted WCA potential, Eq.(6) 
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Parameters. From the phase diagram irf^ we 
set the temperature of our simulations to k^T /£o = 
1.1, which allows for a broad range of (Oc (between 
1.3(7 and 1.7(7) within which the membrane is in 
the fluid state. Furthermore, the bending rigidity 
was calculated as a function of £0c over this range 
of values for kgT /Cq = 1.1 in the same work; re- 
sults from that reference are shown in Figure la. 
Values of the bilayer density calculated in our sim- 
ulations over a similar range of COc are shown in 
Figure lb. 

The units of energy, length, and time in our sim- 
ulations are respectively Co, (7 and Tq. The remain- 
ing parameters can be assigned physical values by 
setting the system to room temperature, T = 300^, 
and noting that the typical width of a lipid bilayer 
is around 5 nm, and the mass of a typical phospho- 
lipid is about 660 g/mol. The units of our system 
can then be assigned as follows: o = 0.9 nm, niQ = 



220 g/mol, Co = 3.77 x lO^^i j = 227gA^/ps2mol, 
and To = Oy^nio/e = 8.86 ps. 

Simulations 

Molecular dynamics (MD) simulations of budding 
were performed at constant temperature and pres- 
sure using the velocity Verlet algorithm, with a 
Langevin thermostar^ to maintain constant tem- 
perature and a modified Andersen barostatl^ to 
maintain constant membrane tension to represent 
wrapping by an infinite membrane. The time step 
was At = 0.01 To^ the friction constant was 7 = 
Xq^, the box friction for the Andersen barostat 
was %ox = 2 ■ lO^'^ and the box mass Q = 10^^ 
in the system units. The reference pressure, Pq, 
is set to 0, to simulate a tensionless membrane. 
The tension equals the pressure because the nor- 
mal component to the membrane, the z-axis in our 
case, is free to fluctuate and does not contribute 
to the pressure. The x and y components of ve- 
locities and positions are rescaled according to the 
changes in the volume. In order to simulate an 
infinite membrane, periodic boundary conditions 
were employed. 

For most simulations the membrane was com- 
prised of n = 21,492 beads. An initial bilayer 
configuration was relaxed by MD and then placed 
normal to the z-axis in a cubic box of side-length 
L = 63.5(7. The particle was introduced in the cen- 
ter of the box with its pole located about 5 a below 
the membrane surface with zero initial velocity. 

Since the membrane was kept tensionless by 
the barostat, the size of the box decreased dur- 
ing simulations as the particle was wrapped. To 
ensure that there were no finite size effects, addi- 
tional sets of simulations were performed, follow- 
ing the same protocol, for membranes with n = 
48,600 beads and initial box size of 100x100x60 
(7^ and with n = 86,400 beads and initial box 
size 130x130x60 o^. Except where mentioned 
otherwise, results are shown for the system with 
n = 21,492 beads. 

Free energy calculations. In addition to per- 
forming dynamical simulations of budding, we 
calculated the potential of mean force as a func- 
tion of particle penetration using umbrella sam- 
pling. Simulations were performed in which the 
system was biased toward particular values of the 
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penetration p by introducing a biasing function 
f/bias({r}) = \Kumh{p{{r}) - PqY . Here p{{r}) 
is the penetration for a configuration {r} and is 
defined as the distance between the top of the par- 
ticle and the center of mass of the membrane. A 
series of windows were performed at different val- 
ues of pq; for all windows JCumb = 200eo/a^. The 
simulations were started for an unwrapped particle 
{p = —6(7), and initial coordinates for each sub- 
sequent window were obtained from simulations 
in the previous one. Statistics from each window 
were stitched together and re- weighted to obtain 
the unbiased free energy using the weighted his- 
togram analysis method. 

Elastic model 

To evaluate the results of the dynamical simula- 
tions and free energy calculations, we compare the 
simulation results to a simplified elastic model for 
invagination of a particle in a membrane. Our elas- 
tic model closely follows that of Ref.^^ but we 
consider an infinite tensionless membrane rather 
than a vesicle. 

The total energy of the particle-membrane sys- 
tem arises from the energy of adhesion between 
the particle and the membrane (ead) and the elas- 
tic energy of the membrane (em)- Following the 
simulation model, we assume that adhesion is me- 
diated by short-range interactions with energy per 
area — £* so that the total energy of adhesion is 
Cad = — £*«wrap with flwrap the area of the mem- 
brane in contact with the particle. Note that e* 
actually describes a free energy since it includes 
the effects of counterion dissociation and other en- 
tropic factors involved in particle associations, but 
following Ref. we refer to it and the elastic 
terms described next as energies to emphasize that 
we are neglecting the (small) contribution to the 
free energy associated with fluctuations around the 
lowest free energy membrane configuration. 

To calculate the elastic contributions to the en- 
ergy, we consider the Helfrich Hamiltonian for an 
infinitesimally thin membrane*^ 

e^ = jda(^Os + ^{2H-Cof + KqK^ (7) 

where cr, is the surface tension, K and Kq are the 
bending rigidity and the Gaussian curvature mod- 



ulus respectively, H and K are the mean and Gaus- 
sian curvatures, and Co is the spontaneous curva- 
ture. Our model membrane is symmetric and ten- 
sionless, so Co and Os are 0. We will use this elas- 
tic model to describe the budding process up un- 
til the point of scission at the neck, and thus the 
topology of the membrane remains constant. As- 
suming that the Gaussian curvature modulus Kg is 
invariant throughout the membrane, the last term 
in Eq.(7) is constant under the Gauss-Bonnet the- 
orem.f^Sl in practice, the properties of the mem- 
brane and thus Kq could change in the vicinity of 
the adsorbed particle, but this effect contributes a 
factor proportional to the adhesion area awrap and 
thus only renormalizes the adhesion free energy 
£*. The elastic energy for a general configuration 
of the membrane is then given by 



K / 1 1 
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with ri and r2 the principal radii of curvature. 




Figure 2: Cross-section of the 3D geometry used for 
the elastic model of a membrane wrapping a particle of 
radius R. The particle, depicted as a green sphere, sticks 
to a section of the membrane in red with area a wrap • 
The surrounding membrane, with area anm, drawn in 
blue, decays toward the flat configuration. The shape 
of this surrounding membrane is taken to be a section 
of a torus for simplicity. 5 and p stand for the outer 
and inner radius of the section of the torus formed by 
the rim region, and a and (j) represent the polar and 
azimuthal angle of spherical coordinates. For a given 
penetration, p, there is a wrapping degree, 6, that min- 
imizes the elastic energy. 

We follow Desemo et al.'^ to assume that the 
following geometry closely corresponds to the 
lowest free energy configuration for a partially 
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wrapped particle (Figure 2). There is an area awrap 
of the membrane tightly adhered to the particle 
with radius of curvature approximately equal to 
the particle radius R, and a rim area, ai-im, between 
the point at which the membrane separates from 
the particle and where it recovers a flat configura- 
tion. Because the particle is a featureless sphere, 
we assume that the lowest free energy configura- 
tion is axisymmetric, to give the energy 
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where a^j^ is the area of the rim surrounding the 
particle, and ri and are the principal radii of cur- 
vature in the rim area. 

We now recast Eq.(9) in terms of two new vari- 
ables, the latitudinal degree of wrapping 9 and 
the penetration p, which is the distance the par- 
ticle travels along the direction normal to the flat 
membrane, measured from the point at which the 
surface of the particle first touches the flat mem- 
brane. Note that the theoretical penetration and the 
one from simulations have different definitions, al- 
though qualitatively both describe the system in a 
similar way. In Figure 2, a schematic of the system 
is depicted as a 2-D cross-sectional cut, in which 
the red line represents the section of the membrane 
bound to the particle and the blue line represents 
the rim region. As shown in the schematic, we 
assume that the rim corresponds to a section of a 
torus (appearing as a circular arc in the 2-D cross- 
section). Although this is only one of the multiple 
shapes the rim can form, it was shown to closely 
correspond to solutions from a full variational cal- 
culation in Ref.'^ and allows us to write the ge- 
ometric properties of the system as explicit func- 
tions of our parameters. In particular, the radius of 
the torus depends uniquely on the particle size R, 
the wrapping degree 9 and the penetration p\ the 
area element on a torus and the two principal radii 
of curvature are^^ Acinm = — ps'ma)dad^, 
ri = p and r2 = ~ ^ gfn^"" ' where a and ^ are 
the polar and azimuthal angle in spherical coor- 
dinates. With the new parametrization, the area 
of the membrane in contact with the particle for 
a wrapping degree 9, turns out to be awrap(^) = 



2;ri?2(i_cos(0)). 

Therefore, the energy of the system can be writ- 
ten in the following way: 
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da 



where A = AnR^ is the surface area of the parti- 
cle. For a given bending rigidity, particle size and 
membrane -particle interaction, the energy only de- 
pends on the penetration p and the wrapping de- 
gree 9. Finally, for each value of p we mini- 
mize the energy Eq.(lO) with respect to 9 to ob- 
tain the membrane configuration and correspond- 
ing energy as a function of penetration alone. The 
results of the minimization are described below in 
section Phase Diagrams 

Results 

System behavior 

To understand the influence of membrane and par- 
ticle properties on budding, we began by perform- 
ing dynamical simulations for a range of particle- 
membrane interaction strengths, e, particle radius 
R, and tt)c, which controls the areal density of 
lipids, the bending rigidity K, and diffusion rates 
within the membrane, as described in section The 
membrane model. Different values of these pa- 
rameters lead to dramatically different behaviors, 
as shown in phase diagrams presented below (Fig- 
ure 12). First we note that the behaviors can be 
grouped into four classes, which we illustrate by 
describing trajectories observed for various val- 
ues of e and constant values of K = 13.9^67" and 
R=\2o (Figure 3). 

For weak adhesion strengths £, no wrapping oc- 
curs; the membrane continues to exhibit only the 
usual spectrum of thermal fluctuations (Figure 4) 
after the particle adheres to it, and the penetration 
oscillates around negative values (Figure 3, case 
for e = O.Sco). 

For a narrow intermediate range of e, the parti- 
cle adheres to the membrane, but wrapping ceases 
at a partially wrapped state (Figure 5), after which 
the degree of particle penetration into the mem- 
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Figure 3: Particle penetration into the membrane, p, 
as a function of time for molecular dynamics trajecto- 
ries with different values of the adhesion strength for 
a system with particle radius R = 12a and membrane 
bending rigidity K = 13.9^b7"- For e = 0.5eo (red line) 
no wrapping occurs. For e = l.Oeo (green line) and 
£ = 1.2£o (orange line), budding becomes stalled at 
a partially wrapped state whose value increases with 
e. For £ = 2.0eo (blue line), 3.0 (pink hne), and 4.0 
(cyan line) the particle undergoes complete encapsula- 
tion. On the right, slices of the system for correspond- 
ing values of p are shown. Images were generated using 
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Figure 4: Adhesion without wrapping. Slices of con- 
figurations extracted from MD simulations with e = 
leo ,R = 10(7, and k = 13.%b7^- Times shown are 
(a) f = 5 • 10^ To, (b) f = 4 • lO^^ToXc) f = 6 • 10"^ To, (d) 
f = MO^To. 



B.OCo), in which the particle is completely encap- 
sulated (Figure 6). In this case, wrapping proceeds 
steadily until the particle is completely surrounded 
by membrane except for a narrow neck region 
(Figure 6d). Wrapping is then completed when 
a thermal fluctuation causes the neck to break 
and have its sides fused (Figure 6e), after which 
the fully wrapped particle diffuses away from the 
membrane (Figure 6f). Since fusion is a stochas- 
tic event, the budding time can be variable and we 
have observed neck configurations lasting between 
500 and 5000 Xq. The elastic theory predicts that 
the shape and length of the neck depend on the bal- 
ance between the adhesion energy and the bending 
energy with strong adhesion favoring a short neck 
and large bending energies favoring a long neck. 
The simulation results are consistent with this pre- 
diction; example configurations are shown in Fig- 
ure 7. The figure shows snapshots from simula- 
tions with particle radius R = 6a, bending rigidity 
K = 13.9^37" and different values of e. A small 
particle size was chosen for the figure because the 
relationship between neck configuration and adhe- 
sion energy is most easily visualized when high 
membrane curvature is required for wrapping. The 
fact that fusion is accessible within the course of 
a typical simulation is an interesting contrast be- 
tween the model studied here and that studied by 
Smith et al.,^^ where fusion was observed only for 
inhomogeneous membranes. 
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brane fluctuates around a steady value (Figure 3, 
case for e = 1 .OCq). The average value of the pene- 
tration remained unchanged for as long as we sim- 
ulated (up to 4- 10"^ To)- The final degree of pene- 
tration in this arrested state increases with e, un- 
til approximately the point at which the particle is 
half wrapped (Figure 3, case for e = 1.2£q). A 
further increase in e results in the next class of 
trajectories (Figure 3, case for e = I.Oeq and e = 



Figure 5: Long-lived partial wrapping. Slices of 
configurations extracted from MD simulations with e = 
1.3£o, R = 10a, and k = UMbT. The particle re- 
mains partially wrapped for the length of the simulation 
{t=4- lO^To). (a) f = 5 • lO^To, (b) t = \.5- 104to,(c) 
? = 2.5 • lO'^To, (d) f = 3 • lO'^To. 

For higher values of e wrapping proceeds ex- 
tremely rapidly (Figure 3, case for e = 4.0eo). 
as there is a strong driving force to increase the 
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(d) 
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Figure 6: Wrapping. Slices of configurations at dif- 
ferent times extracted from MD simulations with e = 
3eo,R = 10a and k = 13.9^6^. The membrane wraps 
the particle (a-d) until a neck or channel connecting the 
flat bilayer and the membrane surrounding the particle 
forms (e). Thermal fluctuations break this narrow neck, 
resulting in the encapsulated particle escaping from the 
membrane (f). Configurations are shown for times (a) 
t = 5- 10^ To, (h) t = I- lO'^ToXc) ? = 1.5 • 10"^ To, (d) 
f = 1.95 • 10"^ To, (e) f = 2 • 104To,(f) t = 2.5- lO'^To- 

number of head-particle interactions (Figure 8). 
As the curvature of the membrane in the vicin- 
ity of the wrapping front increases, the membrane 
structure undergoes ruptures in that region (Fig- 
ure 8d), and a pore forms in the membrane (Fig- 
ure 8e). The fully encapsulated particle then dif- 
fuses away and the pore heals through thermal mo- 
tions of the lipids. The formation of a pore dur- 
ing these budding trajectories resembles the pro- 
cess by which a hydrophobic nanoparticle passes 
th rough membranes in the simulations described 
iUjI^ZESIbut the physical driving forces are different 
in this case and the pore arises for kinetic reasons. 
Namely, the collective wrapping process proceeds 
more slowly than ruptures form in the membrane 
due to the large driving force to increase particle- 
head group contacts. 

Free energy calculations. We were particularly 
interested in the partially wrapped states seen in 
the dynamical simulations described in the previ- 
ous section, as the elastic model predicts only fully 
wrapped or non- wrapped states. To determine 
whether or not these observations corresponded 
to equilibrium configurations, the free energy was 
calculated as a function of the penetration using 
umbrella sampling (section System Model). Cal- 
culated free energy projections are shown for three 
values of e in Figure 9, for which the finite-time 
dynamical simulations respectively ended in no 
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Figure 7: The neck profile depends on the adhe- 
sion strength. Membrane configurations are shown 
shortly before the completion of budding for R = 6(7, 
K = 13.9^b7" and different adhesion strengths, (a) A 
relatively small adhesion strength, e = Seq, leads to a 
long neck, (b) For a = 4£o the neck is shorter (c) 
For £ = 5£o, close to the adhesion strength that leads to 
membrane rupture, the neck length is comparable to the 
height of typical membrane fluctuations. The top row 
of images shows a side view of system configurations 
and the bottom row of images gives the corresponding 
side-view slices. 
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Figure 8: Wrapping via membrane rupture. Slices 
of configurations at different times extracted from MD 
simulations with e = Ssq ,R = 10a, and k = l3.9kBT. 
Particle wrapping (a, b) leads to the formation of a 
pore (c, d). Eventually, the enveloped particle leaves 
the membrane (e) and the pore closes (f). Configu- 
rations shown occured at times (a.) t = 2 ■ lO^To, (b) 
t = 5- 10^To,(c) t = 5.4- lO^To, (d) t = 5.6- lO^ToXe) 
t = S-lO^To, (f) f = 9.5 • lO^To. 
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wrapping (e = O.SCo, Figure 9a), partial wrap- 
ping (e = Co, Figure 9b), and complete wrapping 
(e = 1 .25eo, Figure 9c). For the cases of full wrap- 
ping and no wrapping, the calculated free energy 
projections are consistent with the dynamics re- 
sults. Namely, for £ = O.SEq the minimum free 
energy value corresponds to no wrapping with a 
steep penalty for increasing penetration, while for 
£ = 1.25£o the free energy decreases monotoni- 
cally with increasing penetration until the particle 
is completely wrapped. 
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Figure 9: Free energy profiles as a function of the pen- 
etration, p, calculated from MD simulations using um- 
brella sampling are shown for R = 12(7, K = 13.9^b7^> 
and indicated values of e. At the bottom, slices of the 
system as a function of the penetration are shown. 

In contrast, the minimum value in the free en- 
ergy profile for £ = 1.0£o does not correspond to 
the partially wrapped state observed in the dynam- 
ical simulations, but rather corresponds to com- 
plete wrapping. Further comparison of umbrella 
sampling results to dynamical trajectories suggests 
that, while the penetration coordinate p is capable 
of describing the free energy basins correspond- 
ing to the unwrapped and wrapped states, p alone 
is not sufficient to completely describe the tran- 
sition dynamics. I.e., p is a suitable order pa- 
rameter for determining the free energies of the 
stable states, but not a complete reaction coordi- 
nate.^^ To see this, we chose a set of configura- 
tions from the umbrella sampling trajectories with 
different values of p. For each such configura- 
tion we performed several unbiased MD trajecto- 
ries initialized with velocities using different ran- 
dom number seeds to obtain a crude estimate of 



the commitment probability.'^ We found that tra- 
jectories initiated from configurations with small 
values of p < 5c7 fluctuate around that value and 
configurations with p > 18.8c7 progressed steadily 
to complete wrapping. However, configurations 
with 5o < p < 18. Sc tended to fluctuate around 
the value of p corresponding to their initial config- 
uration, which is inconsistent with the free energy 
profile for p > 15(7 and indicates the presence of 
an additional slow degree of freedom at moderate 
e. 

It is not necessary to identify a perfect reaction 
coordinate to fulfill our primary objective of un- 
derstanding the phase behavior, but we did attempt 
to identify the second relevant dynamical degree of 
freedom. Analysis of umbrella sampling configu- 
rations during the equilibration phase of the cal- 
culation indicates that, when the particle is held 
at a fixed penetration, the membrane configura- 
tion gradually relaxes to a state of increased ad- 
hesion and bending (Figure 10). Thus we expect 
that a reaction coordinate capable of describing the 
transition dynamics needs to include an additional 
collective variable that describes adhesion and/or 
membrane bending. Investigating this possibility, 
however, is beyond the scope of the present work 
focused on the phase behavior. 

Phase diagrams 

Based on the results of dynamical simulations 
over a wide range of parameters, as well as um- 
brella sampling at parameter sets near the tran- 
sition between no wrapping and wrapping, we 
determined phase diagrams as functions of ad- 
hesion energy £, membrane bending modulus 
K, and particle radius R (Figure 12). To en- 
able comparison with the elastic theory Eq.(9), 
it is essential to note the the theoretical param- 
eter £* corresponds to the adhesion free energy 
rather than simply the depth of the head group- 
particle attractive potential well £. Therefore, 
we plot the data as a function of the adhesion 
free energy per area, calculated as £*/k^T ~ 
1 + dr ( V-ticie-head('-) _ J j ^ Flgurc 



-T]log 

1 1 . Here we have neglected the tiny contribution 
from cutting off the potential, we assume that each 
lipid head group in contact with the particle ap- 
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time/lO^r 

Figure 10: Evolution of the membrane area in con- 
tact with the particle (a) and bending energy (b) as 
a function of time for MD simulations with a bias- 
ing potential holding the penetration near p = 16.9a, 
[/bias = 0.5K-umb(^' " 16.875)2, and K„^^ = lOOeo/a^. 
The time courses are averaged over 20 independent tra- 
jectories 

proaches roughly along a radial coordinate, and 
we neglect configurational entropy losses endured 
by the lipid molecules during adhesion. The pa- 
rameter sets for which non-wrapping is the equi- 
librium configuration are shown with * symbols, 
while the parameter sets which lead to equilib- 
rium wrapping are separated into those which in- 
volve long-lived partially wrapped structures (x 
symbols), complete wrapping (-1- symbols), and 
those for which the membrane undergoes rupture 
prior to budding (□ symbols). The wrapping bin- 
odal predicted by the elastic theory is shown as a 
dashed line on each plot. We see that while the the- 
ory and simulations agree to within about 0.2^37". 
the theoretical binodal is below the computational 
results. This discrepancy could occur because we 
have not accounted for the configurational entropy 
contributed by the lipids during adhesion or due to 
the fact that the theory assumes an infinitesimally 
thin membrane. 

The complete phase diagram predicted by the 
elastic theory is shown in Figure 13. Here, the 
dashed line is the binodal, given by £* = below 




Figure 1 1 : Relationship between the adhesion en- 
ergy £ and the adhesion free energy per area e* for 
K= BMbT. 

which wrapping is energetically unfavorable and 
the solid line denotes the spinodal, above which 
wrapping proceeds without any energetic barrier. 
In between the lines there is a barrier to wrapping 
which begins at p = 0, meaning that there are no 
long-lived partially wrapped states consistent with 
this theory for infinite membranes. 

Discussion 

Our dynamical simulations of a minimal molecu- 
lar model for the process of passive endo- or ex- 
ocytosis identified four classes of behaviors re- 
sulting from the interaction of a particle with a 
membrane, no wrapping or adhesion resulting in 
a minimal perturbation of membrane configura- 
tions, partial wrapping, complete wrapping, and 
wrapping via rupture of the membrane. Equilib- 
rium calculations showed that there are only two 
equilibrium configurations, corresponding to no 
wrapping or complete wrapping, and this equi- 
librium phase behavior for the molecular model 
shows strong agreement with the predictions of 
a simplified elastic theory.^ The primary differ- 
ence between the elastic model and the finite-time 
dynamical simulation results is the appearance of 
long-lived partial wrapping states. 

Since the long-lived partially wrapped states 
seen in this study could be significant for dynami- 
cal, time-sensitive particle uptake processes such 
as endocytosis or viral budding in living organ- 
isms, it is worth comparing them to observations 
of other models. Most closely related to our re- 
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Figure 12: Phase diagrams obtained from MD simu- 
lations and free energy calculations (a) as a function of 
£ and K for constant R = 10a and (b) as a function of 
tlie particle size R and adhesion free energy e* for con- 
stant bending rigidity, K = l3.9kBT. Parameter sets are 
identified as those which lead to no wrapping (* sym- 
bols), long-lived partially wrapped structures (x sym- 
bols), complete wrapping (-1- symbols), and those for 
which the membrane undergoes rupture prior to bud- 
ding (□ symbols). The wrapping binodal predicted by 
the elastic theory is shown as a dashed line on each plot. 
The relationship between the adhesion free energy e* 
and the head group-particle attractive well depth e is 
given in the text and in Figure 1 1 . 
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non wrapping 
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Figure 13: Phases diagrams obtained from the energy 
minimization in Eq.(lO) (a) as a function of e* and K 
for constant R = 10(7 and (b) as a function of the par- 
ticle size, R, and the adhesion free energy per area, £*, 
for constant bending rigidity K = 13.%b7"- The bin- 
odal (£* = ^) above which wrapping is energetically 
favorable is indicated by a dashed line and the spinodal, 
above which wrapping proceeds without an energy bar- 
rier, is shown by a solid line. 
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suits, Yue and Zhangl^ study a model compara- 
ble to ours except that the particle is coated with 
ligands that experience attractive interactions with 
membrane lipids. It appears that some of the con- 
figurations which they denote as 'adhesion' corre- 
spond to our long-lived partially wrapped states; 
however, it would be necessary to perform free 
energy calculations to determine whether they are 
metastable as we find here. In contrast to their 
model, we do not identify any parameter sets for 
which the particle partially penetrates into the hy- 
drophobic interior of the membrane. Partially 
wrapped states have also been predicted from 
equilibrium theories in the context of finite-sized 
systems. Deserno et a\J^ examined the budding 
of colloid from the interior of a spherical vesicle 
using the same Hamiltonian as introduced in elas- 
tic theory described above. They found that par- 
tial wrapping corresponds to the equilibrium state 
when the vesicle size is on the order of the col- 
loid diameter due to the increase in curvature en- 
ergy of the finite- sized vesicle. We similarly ob- 
tain partial wrapping configurations as equilibrium 
solutions if we introduce finite size into the elas- 
tic theory studied here (Eq.(9)) by minimizing the 
energy of the elastic theory (Eq.(lO)) with the to- 
tal membrane surface area constrained to 16R^, so 
that the wrapping area plus the rim area cannot 
exceed the total area. Zhang and Nguyen also 
identify partially wrapped states as equilibrium so- 
lutions to an elastic theory, but they observed that 
the catenoid configuration is the only solution to 
the full variational problem for a tensionless infi- 
nite membrane. This solution implies that the elas- 
tic energy of the rim is always zero, and thus wrap- 
ping is only determined by the balance between 
the bending energy in the wrapped region and the 
adhesion energy, which does not lead to partial 
wrapping states. Because the toroid approxima- 
tion for membrane configurations assumed in our 
elastic model is more restrictive than the full vari- 
ational problem, the wrapping binodal shown in 
Figure 13 is shifted to slightly higher values of 
the adhesion energy e than obtained for their the- 
ory,'^ but the behavior is qualitatively unchanged. 
Importantly, neither theory predicts the partially 
wrapped state as a metastable configuration in an 
infinite membrane. 
The membrane size in our simulations was 



chosen large enough to ensure that the theoret- 
ically predicted finite-size effects would not af- 
fect our results. To confirm that this was the 
case, we ran additional dynamical simulations and 
umbrella sampling calculations with membranes 
which were 50% and 100% larger (16200 and 
28800 lipids respectively). The simulation results 
were the same for all three membrane sizes. 

Finally, we consider our minimal model for pas- 
sive endo- or exocytosis in the context of physical 
systems. Based on the length scales discussed in 
section Parameters, the particle diameters in our 
simulations range from 9 to 36 nm. Nanoparti- 
cles are available in a wide range of sizes, with 
particles smaller than 50 nm undergoing the most 
efficient uptake,!^ and our simulation results in- 
dicate a range of possible uptake pathways. Our 
simulated particles are somewhat smaller than the 
size of viral capsids that undergo budding, which 
range from about 40 nm (e.g. hepadnavirusl^ to 
hundreds of nanometers, but the results can be ex- 
trapolated into that range. As shown in Figure 13 
the adhesion-wrapping transition decreases with 
radius as 1/7? for constant bending rigidity. 

Our simulation results indicate that the existence 
of attractive interactions between a particle and 
lipid head groups, which has been propos ed as 
the minimal requirement for viral budding,'^2EI] 
is indeed sufficient to drive efficient wrapping. 
However, to avoid stalled partially-wrapping states 
and membrane rupture, the system would be con- 
fined to a relatively narrow range of adhesion 
strengths spanning about 2k^T / (the -i- sym- 
bols in Figure 12). While this result is qualita- 
tive, since the range increases in width with the 
particle size and the three-bead representation of 
the lipid molecule may lead to model membranes 
which are more susceptible to rupture than those 
comprised of a more realistic lipid, it does estab- 
lish important constraints on viral evolution if bud- 
ding is limited to these ingredients. However, de- 
pending on the viral system, a number of addi- 
tional phenomena contribute to budding, including 
membrane-associated viral envelope or spike pro- 
teins, preferential budding from lipid rafts, the 
use of cytosk eleta l machiner y to actively drive or 
assist buddingl^ or scission,L^'^and the ability of 
the virus to remodel cell membrane properties.'^ 
These effects can broaden the range of functional 
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adhesion energies; e.g. using actin to drive as- 
sembly and budding^ could enable efficient viral 
egress even for adhesion energies at which spon- 
taneous dynamics become stalled. In this case the 
presence of a barrier to budding could serve as a 
regulatory feature. The agreement between our 
simulation results and the elastic theory over some 
ranges of parameter space indicates that some or 
all of these effects could in principle be captured 
by extending existing elastic theories along the 
lines of Liu et. al.'s description of active endo- 
cytosiSjI^but care would be required to include all 
relevant slow degrees of freedom near transitions 
between wrapping and no wrapping. 
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